{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "execution": {},
    "heading_collapsed": false,
    "level": 1
   },
   "source": [
    "# Biological Neuronal Networks: The Leaky Integrate-and-Fire (LIF) Neuronal Network\n",
    "\n",
    "__Referred from:__ [Brian2 documantations](https://brian2.readthedocs.io/en/stable/)\n",
    "\n",
    "__Content modified:__ Kai Chen \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "execution": {}
   },
   "outputs": [],
   "source": [
    "#@title load matplotlib and rcParams\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "nma_style = {\n",
    "    'figure.figsize' : (8, 6),\n",
    "    'figure.autolayout' : True,\n",
    "    'font.size' : 15,\n",
    "    'xtick.labelsize' : 'small',\n",
    "    'ytick.labelsize' : 'small',\n",
    "    'legend.fontsize' : 'small',\n",
    "    'axes.spines.top' : False,\n",
    "    'axes.spines.right' : False,\n",
    "    'xtick.major.size' : 5,\n",
    "    'ytick.major.size' : 5,\n",
    "}\n",
    "for key, value in nma_style.items():\n",
    "    plt.rcParams[key] = value\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "execution": {}
   },
   "source": [
    "## Import functions and classes in Brian2\n",
    "\n",
    "All Brian scripts start with the following."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {}
   },
   "outputs": [],
   "source": [
    "from brian2 import *"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "execution": {},
    "heading_collapsed": false,
    "level": 2
   },
   "source": [
    "## Units system\n",
    "\n",
    "Brian has a system for using quantities with physical dimensions. \n",
    "\n",
    "The base units are defined by their standard **SI** unit names: `amp`/`ampere`, `kilogram`/`kilogramme`, `second`, `metre`/`meter`, `mole`/`mol`, `kelvin`, and `candela`. \n",
    "\n",
    "In addition to these base units, Brian defines a set of derived units: `coulomb`, `farad`, `gram`/`gramme`, `hertz`, `joule`, `liter`/`litre`, `molar`, `pascal`, `ohm`, `siemens`, `volt`, `watt`, together with prefixed versions (e.g. `msiemens` = `0.001*siemens`) using the prefixes `p`, `n`, `u`, `m`, `k`, `M`, `G`, `T` (two exceptions to this rule: `kilogram` is not defined with any additional prefixes, and `metre` and meter are additionaly defined with the “`centi`” prefix, i.e. `cmetre`/`cmeter`). \n",
    "\n",
    "For convenience, a couple of additional useful standard abbreviations such as `cm` (instead of `cmetre`/`cmeter`), `nS` (instead of `nsiemens`), `ms` (instead of `msecond`), `Hz` (instead of `hertz`), `mM` (instead of `mmolar`) are included. To avoid clashes with common variable names, no one-letter abbreviations are provided (e.g. you can use `mV` or `nS`, but not `V` or `S`).\n",
    "\n",
    "Let's try a few of them!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {}
   },
   "outputs": [],
   "source": [
    "20*volt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {}
   },
   "outputs": [],
   "source": [
    "1000*amp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {}
   },
   "outputs": [],
   "source": [
    "1e6*volt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {}
   },
   "outputs": [],
   "source": [
    "1000*namp"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "execution": {}
   },
   "source": [
    "Also note that combinations of units with work as expected:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {}
   },
   "outputs": [],
   "source": [
    "10*nA*5*Mohm"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "execution": {}
   },
   "source": [
    "And if you try to do something wrong like adding amps and volts, what happens?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": false,
    "level": 2
   },
   "source": [
    "## Build A Leaky Integrate-and-Fire Neuron from the ground\n",
    "\n",
    "In Brian, all models are defined by systems of differential equations. Here's the first simple example:\n",
    "\n",
    "consider the LIF neuron defined as:\n",
    "\n",
    "$$\n",
    "C_m\\frac{dv}{dt} = -g_L(v-E_L) + I_{ext}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {}
   },
   "outputs": [],
   "source": [
    "start_scope()\n",
    "\n",
    "# define neuronal parameters\n",
    "C_m = 0.5*nF\n",
    "g_L = 50*nS\n",
    "E_L = -75*mV\n",
    "V_rest = E_L\n",
    "I_ext = 1*nA\n",
    "\n",
    "# using multi-line string to define the system of ODE for our LIF model\n",
    "eqs = '''\n",
    "dv/dt = (-g_L*(v-E_L) + I_ext)/ C_m : volt\n",
    "'''\n",
    "\n",
    "G = NeuronGroup(1, eqs, method='exact')\n",
    "G.v = V_rest    # set initial value of v\n",
    "run(30*ms)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "execution": {}
   },
   "source": [
    "In Python, the notation ``'''`` is used to begin and end a multi-line string. So the equations are just a string with one line per equation. The equations are formatted with standard mathematical notation, with one addition. At the end of a line you write ``: unit`` where ``unit`` is the SI unit of that variable.\n",
    "Note that this is not the unit of the two sides of the equation (which would be ``1/second``), but the unit of the *variable* defined by the equation, i.e. in this case $v$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "execution": {},
    "level": 7
   },
   "source": [
    "First off, ignore that ``start_scope()`` at the top of the cell. You'll see that in each cell in this tutorial where we run a simulation. All it does is make sure that any Brian objects created before the function is called aren't included in the next run of the simulation.\n",
    "\n",
    "Secondly, you'll see that there is an \"INFO\" message about not specifying the numerical integration method. This is harmless and just to let you know what method we chose, but we'll fix it in the next cell by specifying the method explicitly.\n",
    "\n",
    "So, what has happened here? Well, the command ``run(30*ms)`` runs the simulation for 30 ms.\n",
    "\n",
    "Now let's take a look at a graph of how the variable ``v`` evolves over time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {}
   },
   "outputs": [],
   "source": [
    "start_scope()\n",
    "\n",
    "# defineriable monitor to record state of v\n",
    "G = NeuronGroup(1, eqs, method='exact')\n",
    "G.v = V_rest\n",
    "M = StateMonitor(G, 'v', record=True)\n",
    "\n",
    "print('Before simulation: v = %5.2f mV' % (G.v[0]/mV))\n",
    "run(100*ms)\n",
    "print('After simulation:  v = %5.2f mV' % (G.v[0]/mV))\n",
    "\n",
    "# plot the trace of V\n",
    "plt.plot(M.t/ms, M.v[0]/mV)\n",
    "plt.xlabel('Time (ms)')\n",
    "plt.ylabel('v (mV)');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "execution": {}
   },
   "source": [
    "In this example, we used the object ``StateMonitor`` object. This is used to record the values of a neuron variable while the simulation runs. The first two arguments are the group to record from, and the variable you want to record from. We also specify ``record=0``. This means that we record all values for neuron 0. We have to specify which neurons we want to record because in large simulations with many neurons it usually uses up too much RAM to record the values of all neurons.\n",
    "\n",
    "Now try modifying the equations and parameters and see what happens in the cell below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {}
   },
   "outputs": [],
   "source": [
    "start_scope()\n",
    "\n",
    "# change external constant drive to sine wave.\n",
    "eqs_sine = '''\n",
    "dv/dt = ( -g_L*(v-E_L) + I_ext*sin(2*pi*100*Hz*t) ) / C_m : volt\n",
    "'''\n",
    "\n",
    "# Change to Euler method because exact integrator doesn't work here\n",
    "G = NeuronGroup(1, eqs_sine, method='euler')\n",
    "M = StateMonitor(G, 'v', record=0)\n",
    "\n",
    "G.v = V_rest # initial value\n",
    "\n",
    "run(60*ms)\n",
    "\n",
    "plt.plot(M.t/ms, M.v[0]/mV)\n",
    "plt.xlabel('Time (ms)')\n",
    "plt.ylabel('v (mV)');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "execution": {},
    "heading_collapsed": false,
    "level": 2
   },
   "source": [
    "## Adding spikes\n",
    "\n",
    "So far we haven't done anything neuronal, just played around with differential equations. Now let's start adding spiking behaviour."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {}
   },
   "outputs": [],
   "source": [
    "start_scope()\n",
    "\n",
    "I_ext = 2*nA    # larger external drive to ensure firing\n",
    "V_th = -55*mV\n",
    "V_reset = V_rest\n",
    "\n",
    "G = NeuronGroup(1, eqs, threshold='v>V_th', reset='v = V_reset', method='exact')\n",
    "G.v = V_rest\n",
    "\n",
    "M = StateMonitor(G, 'v', record=0)\n",
    "spikemon = SpikeMonitor(G)\n",
    "\n",
    "run(50*ms)\n",
    "plt.plot(M.t/ms, M.v[0]/mV)\n",
    "plt.xlabel('Time (ms)')\n",
    "plt.ylabel('v (mV)')\n",
    "\n",
    "print('Spike times: %s' % spikemon.t[:])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "execution": {}
   },
   "source": [
    "We've added two new keywords to the ``NeuronGroup`` declaration: ``threshold='v>V_th'`` and ``reset='v = V_reset'``. What this means is that when ``v>V_th`` we fire a spike, and immediately reset ``v = V_reset`` after the spike. We can put any expression and series of statements as these strings.\n",
    "\n",
    "As you can see, at the beginning the behaviour is the same as before until ``v`` crosses the threshold ``v>V_th`` at which point you see it reset to 0. You can't see it in this figure, but internally Brian has registered this event as a spike. Let's have a look at that."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "execution": {}
   },
   "source": [
    "The ``SpikeMonitor`` object takes the group whose spikes you want to record as its argument and stores the spike times in the variable ``t``. Let's plot those spikes on top of the other figure to see that it's getting it right."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "execution": {}
   },
   "source": [
    "Now try changing the strings for ``threshold`` and ``reset`` in the cell above to see what happens."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "execution": {},
    "heading_collapsed": false,
    "level": 2
   },
   "source": [
    "## Refractoriness\n",
    "\n",
    "A common feature of neuron models is refractoriness. This means that after the neuron fires a spike it becomes refractory for a certain duration and cannot fire another spike until this period is over. Here's how we do that in Brian."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {}
   },
   "outputs": [],
   "source": [
    "start_scope()\n",
    "\n",
    "I_ext = 2*nA    # larger external drive to ensure firing\n",
    "V_th = -55*mV\n",
    "V_reset = V_rest\n",
    "\n",
    "eqs = '''\n",
    "dv/dt = (-g_L*(v-E_L) + I_ext)/ C_m : volt (unless refractory)\n",
    "'''\n",
    "\n",
    "G = NeuronGroup(1, eqs, threshold='v>V_th', reset='v = V_reset', refractory=5*ms, method='exact')\n",
    "G.v = V_rest\n",
    "\n",
    "M = StateMonitor(G, 'v', record=0)\n",
    "spikemon = SpikeMonitor(G)\n",
    "\n",
    "run(50*ms)\n",
    "plt.plot(M.t/ms, M.v[0]/mV)\n",
    "for t in spikemon.t:\n",
    "    plt.axvline(t/ms, ls='--', c='C1', lw=3)   # using matplotlib.pyplot.axvline to indicate spike times\n",
    "plt.xlabel('Time (ms)')\n",
    "plt.ylabel('v (mV)');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "execution": {}
   },
   "source": [
    "As you can see in this figure, after the first spike, ``v`` stays at 0 for around 5 ms before it resumes its normal behaviour. To do this, we've done two things. Firstly, we've added the keyword ``refractory=5*ms`` to the ``NeuronGroup`` declaration. On its own, this only means that the neuron cannot spike in this period (see below), but doesn't change how ``v`` behaves. In order to make ``v`` stay constant during the refractory period, we have to add ``(unless refractory)`` to the end of the definition of ``v`` in the differential equations. What this means is that the differential equation determines the behaviour of ``v`` unless it's refractory in which case it is switched off."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "execution": {},
    "heading_collapsed": false,
    "level": 2
   },
   "source": [
    "## Multiple neurons\n",
    "\n",
    "So far we've only been working with a single neuron. Let's do something interesting with multiple neurons."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {}
   },
   "outputs": [],
   "source": [
    "start_scope()\n",
    "\n",
    "N = 100\n",
    "\n",
    "G = NeuronGroup(N, eqs, threshold='v>V_th', reset='v=V_reset', refractory=2*ms, method='exact')\n",
    "G.v = 'rand()*(V_th-V_rest)+V_rest'     # randomly assign the initial v;\n",
    "\n",
    "spikemon = SpikeMonitor(G)\n",
    "\n",
    "run(50*ms)\n",
    "\n",
    "plt.plot(spikemon.t/ms, spikemon.i, '.k')\n",
    "plt.xlabel('Time (ms)')\n",
    "plt.ylabel('Neuron index');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "execution": {}
   },
   "source": [
    "This shows a few changes. Firstly, we've got a new variable ``N`` determining the number of neurons. Secondly, we added the statement ``G.v = 'rand()'`` before the run. What this does is initialise each neuron with a different uniform random value between 0 and 1. We've done this just so each neuron will do something a bit different. The other big change is how we plot the data in the end.\n",
    "\n",
    "As well as the variable ``spikemon.t`` with the times of all the spikes, we've also used the variable ``spikemon.i`` which gives the corresponding neuron index for each spike, and plotted a single black dot with time on the x-axis and neuron index on the y-value. This is the standard \"raster plot\" used in neuroscience."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "execution": {}
   },
   "source": [
    "## Add Synapses\n",
    "\n",
    "Once you have some neurons, the next step is to connect them up via synapses. We'll explore the simplest possible type of synapse that causes an instantaneous change in a variable after a spike."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {}
   },
   "outputs": [],
   "source": [
    "start_scope()\n",
    "\n",
    "eqs = '''\n",
    "dv/dt = (-g_L*(v-E_L) + I)/ C : volt\n",
    "I : amp\n",
    "C : farad\n",
    "'''\n",
    "\n",
    "G = NeuronGroup(2, eqs, threshold='v>V_th', reset='v = V_reset', method='exact')\n",
    "G.v = V_rest\n",
    "G.I = [2, 0]*nA\n",
    "G.C = [1, 10]*nF\n",
    "\n",
    "# Comment these two lines out to see what happens without Synapses\n",
    "S = Synapses(G, G, on_pre='v += 4*mV')\n",
    "S.connect(i=0, j=1)\n",
    "\n",
    "M = StateMonitor(G, 'v', record=True)\n",
    "\n",
    "run(100*ms)\n",
    "\n",
    "plt.plot(M.t/ms, M.v[0]/mV, label='Neuron 0')\n",
    "plt.plot(M.t/ms, M.v[1]/mV, label='Neuron 1')\n",
    "plt.xlabel('Time (ms)')\n",
    "plt.ylabel('v (mV)')\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "execution": {}
   },
   "source": [
    "There are a few things going on here. First of all, let's recap what is going on with the ``NeuronGroup``. We've created two neurons, each of which has the same differential equation but different values for parameters $I$ and $C_m$. Neuron 0 has ``I=2*nA`` and ``C_m=1*nF`` which means that is driven to repeatedly spike at a fairly high rate. Neuron 1 has ``I=0`` and ``C_m=10*nF`` which means that on its own - without the synapses - it won't spike at all (the driving current I is 0). You can prove this to yourself by commenting out the two lines that define the synapse.\n",
    "\n",
    "Next we define the synapses: ``Synapses(source, target, ...)`` means that we are defining a synaptic model that goes from ``source`` to ``target``. In this case, the source and target are both the same, the group ``G``. The syntax ``on_pre='v += 2*mV'`` means that when a spike occurs in the presynaptic neuron (hence ``on_pre``) it causes an instantaneous change on postsynaptic neuron to happen ``v += 2*mV``. So in total, what this model says is that whenever two neurons in G are connected by a synapse, when the source neuron fires a spike the target neuron will have its value of ``v`` increased by 2*mV.\n",
    "\n",
    "However, at this point we have only defined the synapse model, we haven't actually created any synapses. The next line ``S.connect(i=0, j=1)`` creates a synapse from neuron 0 to neuron 1."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "execution": {}
   },
   "source": [
    "## Adding a weight\n",
    "\n",
    "In the previous section, we hard coded the weight of the synapse to be the value `2 mV`, but often we would to allow this to be different for different synapses. We do that by introducing synapse equations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {}
   },
   "outputs": [],
   "source": [
    "start_scope()\n",
    "\n",
    "G = NeuronGroup(3, eqs, threshold='v>V_th', reset='v = V_reset', method='exact')\n",
    "G.v = V_rest\n",
    "G.I = [2, 0, 0]*nA\n",
    "G.C = [1, 10, 10]*nF\n",
    "\n",
    "# Comment these two lines out to see what happens without Synapses\n",
    "S = Synapses(G, G, 'w : volt', on_pre='v += w')\n",
    "S.connect(i=0, j=[1, 2])\n",
    "S.w = 'j*2*mV'\n",
    "\n",
    "M = StateMonitor(G, 'v', record=True)\n",
    "\n",
    "run(100*ms)\n",
    "\n",
    "plt.plot(M.t/ms, M.v[0]/mV, label='Neuron 0')\n",
    "plt.plot(M.t/ms, M.v[1]/mV, label='Neuron 1')\n",
    "plt.plot(M.t/ms, M.v[2]/mV, label='Neuron 2')\n",
    "plt.xlabel('Time (ms)')\n",
    "plt.ylabel('v (mV)')\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "execution": {}
   },
   "source": [
    "This example behaves very similarly to the previous example, but now there's a synaptic weight variable ``w``. The string ``'w : volt'`` is an equation string, precisely the same as for neurons, that defines a single dimensionless parameter ``w``. We changed the behaviour on a spike to ``on_pre='v_post += w'`` now, so that each synapse can behave differently depending on the value of ``w``. To illustrate this, we've made a third neuron which behaves precisely the same as the second neuron, and connected neuron 0 to both neurons 1 and 2. We've also set the weights via ``S.w = 'j*2*mV'``. When ``i`` and ``j`` occur in the context of synapses, ``i`` refers to the source neuron index, and ``j`` to the target neuron index. So this will give a synaptic connection from 0 to 1 with weight ``2*mV=1*2*mV`` and from 0 to 2 with weight ``4*mV=2*2*mV``."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "execution": {}
   },
   "source": [
    "## Introducing a delay\n",
    "\n",
    "So far, the synapses have been instantaneous, but we can also make them act with a certain delay."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {}
   },
   "outputs": [],
   "source": [
    "start_scope()\n",
    "\n",
    "G = NeuronGroup(3, eqs, threshold='v>V_th', reset='v = V_reset', method='exact')\n",
    "G.v = V_rest\n",
    "G.I = [2, 0, 0]*nA\n",
    "G.C = [1, 10, 10]*nF\n",
    "\n",
    "# Comment these two lines out to see what happens without Synapses\n",
    "S = Synapses(G, G, 'w : volt', on_pre='v += w')\n",
    "S.connect(i=0, j=[1, 2])\n",
    "S.w = 'j*2*mV'\n",
    "S.delay = 'j*2*ms'\n",
    "\n",
    "M = StateMonitor(G, 'v', record=True)\n",
    "\n",
    "run(100*ms)\n",
    "\n",
    "plt.plot(M.t/ms, M.v[0]/mV, label='Neuron 0')\n",
    "plt.plot(M.t/ms, M.v[1]/mV, label='Neuron 1')\n",
    "plt.plot(M.t/ms, M.v[2]/mV, label='Neuron 2')\n",
    "plt.xlabel('Time (ms)')\n",
    "plt.ylabel('v (mV)')\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "execution": {}
   },
   "source": [
    "As you can see, that's as simple as adding a line ``S.delay = 'j*2*ms'`` so that the synapse from 0 to 1 has a delay of 2 ms, and from 0 to 2 has a delay of 4 ms"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "execution": {}
   },
   "source": [
    "## More complex connectivity\n",
    "\n",
    "So far, we specified the synaptic connectivity explicitly, but for larger networks this isn't usually possible. For that, we usually want to specify some condition."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {}
   },
   "outputs": [],
   "source": [
    "start_scope()\n",
    "\n",
    "N = 10\n",
    "G = NeuronGroup(N, 'v:1')\n",
    "S = Synapses(G, G)\n",
    "S.connect(condition='i!=j', p=0.2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "execution": {}
   },
   "source": [
    "Here we've created a dummy neuron group of N neurons and a dummy synapses model that doens't actually do anything just to demonstrate the connectivity. The line ``S.connect(condition='i!=j', p=0.2)`` will connect all pairs of neurons ``i`` and ``j`` with probability 0.2 as long as the condition ``i!=j`` holds. So, how can we see that connectivity? Here's a little function that will let us visualise it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "execution": {}
   },
   "outputs": [],
   "source": [
    "#@title Utility function: visualize_connectivity()\n",
    "def visualise_connectivity(S):\n",
    "    Ns = len(S.source)\n",
    "    Nt = len(S.target)\n",
    "    fig, ax = plt.subplots(1,2, figsize=(10, 4))\n",
    "    ax[0].plot(np.zeros(Ns), np.arange(Ns), 'ok', ms=10)\n",
    "    ax[0].plot(np.ones(Nt), np.arange(Nt), 'ok', ms=10)\n",
    "    for i, j in zip(S.i, S.j):\n",
    "        ax[0].plot([0, 1], [i, j], '-k')\n",
    "    ax[0].set_xticks([0, 1],)\n",
    "    ax[0].set_xticklabels(['Source', 'Target'])\n",
    "    ax[0].set_ylabel('Neuron index')\n",
    "    ax[0].set_xlim(-0.1, 1.1)\n",
    "    ax[0].set_ylim(-1, max(Ns, Nt))\n",
    "    S_mat = np.ones((Nt, Ns))\n",
    "    S_mat[S.j, S.i] = 0\n",
    "    ax[1].pcolor(S_mat, cmap='gray', edgecolor='k')\n",
    "    ax[1].invert_yaxis()\n",
    "    ax[1].set_xlabel('Source neuron index')\n",
    "    ax[1].set_ylabel('Target neuron index')\n",
    "    ax[1].axis('scaled')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {}
   },
   "outputs": [],
   "source": [
    "visualise_connectivity(S)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "execution": {}
   },
   "source": [
    "Please check [here](https://brian2.readthedocs.io/en/stable/resources/tutorials/2-intro-to-brian-synapses.html#more-complex-connectivity) for more complex connectivity."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "execution": {}
   },
   "source": [
    "## Simulate randomly connected network with 100 LIF neurons\n",
    "\n",
    "Now, let's put all those pieces together and simulated a random network with 100 LIF neurons."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {}
   },
   "outputs": [],
   "source": [
    "start_scope()\n",
    "\n",
    "# define neuronal parameters\n",
    "C_m = 0.5*nF\n",
    "g_L = 25*nS\n",
    "E_L = 10*mV\n",
    "V_rest = E_L\n",
    "V_reset = V_rest\n",
    "V_th = 30*mV\n",
    "I_ext = 1*nA\n",
    "\n",
    "# define the system of LIF model\n",
    "eqs = '''\n",
    "dv/dt = (-g_L*(v-E_L) + I_ext)/C_m : volt (unless refractory)\n",
    "'''\n",
    "\n",
    "N = 100\n",
    "\n",
    "G = NeuronGroup(N, eqs, threshold='v>V_th', reset='v=V_reset', refractory=2*ms, method='euler')\n",
    "G.v = 'rand()*(V_th-V_rest)+V_rest'\n",
    "# G.v = V_rest\n",
    "\n",
    "S = Synapses(G, G, on_pre='v+=0.5*mV', delay=2*ms)\n",
    "S.connect(condition='i!=j', p=0.2)\n",
    "\n",
    "spikemon = SpikeMonitor(G)\n",
    "\n",
    "run(50*ms)\n",
    "\n",
    "plt.plot(spikemon.t/ms, spikemon.i, '.k')\n",
    "plt.ylim(0,N)\n",
    "plt.xlabel('Time (ms)')\n",
    "plt.ylabel('Neuron index');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "execution": {}
   },
   "source": [
    "Note that we also changed the `method` keyword argument to use `'euler'` (which stands for the [Euler-Maruyama method](https://en.wikipedia.org/wiki/Euler%E2%80%93Maruyama_method)); the 'exact' method that we used earlier is not applicable to those stochastic differential equations."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "execution": {}
   },
   "source": [
    "\n",
    "## Coding Exercise 1: Simulate LIF network containing both excitatory and inhibitory neurons.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {}
   },
   "outputs": [],
   "source": [
    "start_scope()\n",
    "\n",
    "# USE these parameters. (Of course, you can try other values latter.)\n",
    "C_m = 0.5*nF\n",
    "g_L = 25*nS\n",
    "E_L = 10*mV\n",
    "V_rest = E_L\n",
    "V_reset = V_rest\n",
    "V_th = 30*mV\n",
    "I_ext = 1*nA\n",
    "\n",
    "# define the system of LIF model\n",
    "eqs = '''\n",
    "dv/dt = (-g_L*(v-E_L) + I_ext)/C_m : volt (unless refractory)\n",
    "'''\n",
    "\n",
    "N_E = 400\n",
    "N_I = 100\n",
    "\n",
    "G_E = NeuronGroup(N_E, eqs, threshold='v>V_th', reset='v=V_reset', refractory=2*ms, method='euler')\n",
    "G_I = NeuronGroup(N_I, eqs, threshold='v>V_th', reset='v=V_reset', refractory=2*ms, method='euler')\n",
    "G_E.v = 'rand()*(V_th-V_rest)+V_rest'\n",
    "G_I.v = 'rand()*(V_th-V_rest)+V_rest'\n",
    "\n",
    "S_EE = Synapses(G_E, G_E, on_pre='v+=0.5*mV',   delay=2*ms)\n",
    "S_IE = Synapses(G_E, G_I, on_pre='v+=0.5*mV',   delay=2*ms)\n",
    "S_EI = Synapses(G_I, G_E, on_pre='v-=0.5*mV*4', delay=2*ms)\n",
    "S_II = Synapses(G_I, G_I, on_pre='v-=0.5*mV*4', delay=2*ms)\n",
    "for S in (S_EE, S_IE, S_EI, S_II):\n",
    "    S.connect(condition='i!=j', p=0.1)\n",
    "\n",
    "spikemon_E = SpikeMonitor(G_E)\n",
    "spikemon_I = SpikeMonitor(G_I)\n",
    "\n",
    "run(50*ms)\n",
    "\n",
    "plt.plot(spikemon_E.t/ms, spikemon_E.i, '.b')\n",
    "plt.plot(spikemon_I.t/ms, spikemon_I.i+N_E, '.r')\n",
    "plt.ylim(0,N_E+N_I)\n",
    "plt.xlabel('Time (ms)')\n",
    "plt.ylabel('Neuron index');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {}
   },
   "outputs": [],
   "source": [
    "run(50*ms)\n",
    "\n",
    "plt.plot(spikemon_E.t/ms, spikemon_E.i, '.b')\n",
    "plt.plot(spikemon_I.t/ms, spikemon_I.i+N_E, '.r')\n",
    "plt.ylim(0,N_E+N_I)\n",
    "plt.xlabel('Time (ms)')\n",
    "plt.ylabel('Neuron index');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "execution": {}
   },
   "source": [
    "## Explore different dynamical states of E-I LIF networks\n",
    "\n",
    "**Reference:**\n",
    "Brunel, N. (2000). Phase diagrams of sparsely connected networks of excitatory and inhibitory spiking neurons. Neurocomputing, 32–33, 307–312. https://doi.org/10.1016/S0925-2312(00)00179-X\n",
    "\n",
    "- **Synchronous regular (SR) states**: where neurons are almost fully synchronized in a few clusters and behave as oscillators when excitation dominates inhibition and synaptic time distributions are sharply peaked;\n",
    "- **Asynchronous regular (AR) states**: with stationary global activity and quasiregular individual neuron firing when excitation dominates inhibition and synaptic time distributions are broadly peaked;\n",
    "- **Asynchronous irregular (AI) states**: with stationary global activity but strongly irregular individual firing at low rates when inhibition dominates excitation in an intermediate range of external frequencies; \n",
    "- **Synchronous irregular (SI) states**: with oscillatory global activity but strongly irregular individual firing at low (compared to the global oscillation frequency) firing rates, when inhibition dominates excitation and either low external frequencies (slow oscillations) or high external frequencies (fast oscillations). When the average synaptic time constant is high enough, these two regions merge together."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "execution": {}
   },
   "source": [
    "## Coding Exercise: Tuning parameters to achieve SI state."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {}
   },
   "outputs": [],
   "source": [
    "start_scope()\n",
    "\n",
    "# USE these parameters. (Of course, you can try other values latter)\n",
    "tau_m = 20*ms       # membrane time constant\n",
    "V_reset = 10*mV     # reset potential\n",
    "V_th = 20*mV        # threshold potential\n",
    "J = 0.2*mV          # EPSP\n",
    "g = 6               # ratio of IPSP/EPSP\n",
    "delay=1.8*ms        # synaptic transimition delay\n",
    "ext_ratio = 0.9     # ratio of nu_ext/nu_th\n",
    "p = 0.1             # connectivity probability\n",
    "\n",
    "# define the system of LIF model\n",
    "eqs = '''\n",
    "dv/dt = -v/tau_m : volt (unless refractory)\n",
    "'''\n",
    "\n",
    "ratio = 4\n",
    "N_I = 2500\n",
    "N_E = N_I*ratio\n",
    "\n",
    "# write your own code to construct the network\n",
    "\n",
    "\n",
    "# define external Poisson drive\n",
    "nu_th = V_th/(p*N_E*J*tau_m)\n",
    "Poisson_E = PoissonGroup(N_E, rates=nu_th*N_E*p*ext_ratio)\n",
    "Poisson_I = PoissonGroup(N_I, rates=nu_th*N_E*p*ext_ratio)\n",
    "\n",
    "# Also use Synapses to inject outputs from PoissonGroup to NeuronGroup\n",
    "\n",
    "# write your own code to run simulation and plot the results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {}
   },
   "outputs": [],
   "source": [
    "# to_remove solution\n",
    "\n",
    "start_scope()\n",
    "\n",
    "# USE these parameters. (Of course, you can try other values latter)\n",
    "tau_m = 20*ms       # membrane time constant\n",
    "V_reset = 10*mV     # reset potential\n",
    "V_th = 20*mV        # threshold potential\n",
    "J = 0.2*mV          # EPSP\n",
    "g = 6               # ratio of IPSP/EPSP\n",
    "delay=1.8*ms        # synaptic transimition delay\n",
    "ext_ratio = 0.9     # ratio of nu_ext/nu_th\n",
    "p = 0.1             # connectivity probability\n",
    "\n",
    "# define the system of LIF model\n",
    "eqs = '''\n",
    "dv/dt = -v/tau_m : volt (unless refractory)\n",
    "'''\n",
    "\n",
    "ratio = 4\n",
    "N_I = 2500\n",
    "N_E = N_I*ratio\n",
    "\n",
    "# write your own code to construct the network\n",
    "G_E = NeuronGroup(N_E, eqs, threshold='v>V_th', reset='v=V_reset', refractory=2*ms, method='euler')\n",
    "G_I = NeuronGroup(N_I, eqs, threshold='v>V_th', reset='v=V_reset', refractory=2*ms, method='euler')\n",
    "G_E.v = 'rand()*V_th'\n",
    "G_I.v = 'rand()*V_th'\n",
    "\n",
    "S_EE = Synapses(G_E, G_E, 'w:volt',  on_pre='v+=w', delay=delay)\n",
    "S_IE = Synapses(G_E, G_I, 'w:volt',  on_pre='v+=w', delay=delay)\n",
    "S_EI = Synapses(G_I, G_E, 'w:volt',  on_pre='v-=w', delay=delay)\n",
    "S_II = Synapses(G_I, G_I, 'w:volt',  on_pre='v-=w', delay=delay)\n",
    "\n",
    "\n",
    "for S in (S_EE, S_II):\n",
    "    S.connect(condition='i!=j', p=p)\n",
    "for S in (S_IE, S_EI):\n",
    "    S.connect(p=p)\n",
    "\n",
    "for S in (S_EE, S_IE):\n",
    "    S.w = J\n",
    "for S in (S_EI, S_II):\n",
    "    S.w = J*g\n",
    "\n",
    "# define external Poisson drive\n",
    "nu_th = V_th/(p*N_E*J*tau_m)\n",
    "Poisson_E = PoissonGroup(N_E, rates=nu_th*N_E*p*ext_ratio)\n",
    "Poisson_I = PoissonGroup(N_I, rates=nu_th*N_E*p*ext_ratio)\n",
    "\n",
    "# Also use Synapses to inject outputs from PoissonGroup to NeuronGroup\n",
    "SP_E = Synapses(Poisson_E, G_E, on_pre='v+=J')\n",
    "SP_I = Synapses(Poisson_I, G_I, on_pre='v+=J')\n",
    "SP_E.connect(condition='i==j')\n",
    "SP_I.connect(condition='i==j')\n",
    "\n",
    "# write your own code to run simulation and plot the results\n",
    "\n",
    "# define Monitors\n",
    "spikemon_E = SpikeMonitor(G_E)\n",
    "spikemon_I = SpikeMonitor(G_I)\n",
    "pop_E = PopulationRateMonitor(G_E)\n",
    "pop_I = PopulationRateMonitor(G_I)\n",
    "\n",
    "# run simulation\n",
    "T=500*ms\n",
    "run(T)\n",
    "\n",
    "fig, ax = plt.subplots(2,1, figsize=(15,10))\n",
    "ax[0].plot(spikemon_E.t/ms, spikemon_E.i, '.b', ms=2)\n",
    "ax[0].plot(spikemon_I.t/ms, spikemon_I.i+N_E, '.r', ms=2)\n",
    "ax[0].set_xlim(0, T/ms)\n",
    "ax[0].set_ylim(0,N_E+N_I)\n",
    "ax[0].set_xlabel('Time (ms)')\n",
    "ax[0].set_ylabel('Neuron index')\n",
    "ax[0].set_title(f'Individual firing rate : {(spikemon_E.i.shape[0] + spikemon_I.i.shape[0])/(N_E+N_I)/T/Hz:5.2f} Hz')\n",
    "\n",
    "ax[1].plot(pop_E.t / ms, pop_E.smooth_rate(width=2 * ms) / Hz, color='b', label='excitatory')\n",
    "ax[1].plot(pop_I.t / ms, pop_I.smooth_rate(width=2 * ms) / Hz, color='r', label='inhibitory')\n",
    "ax[1].legend()\n",
    "ax[1].set_xlim(0, T/ms)\n",
    "ax[1].set_xlabel('Time (ms)')\n",
    "ax[1].set_ylabel('Populational Firing Rate (Hz)');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "execution": {},
    "heading_collapsed": false,
    "level": 2
   },
   "source": [
    "Once you got the network in **SI** state, you may also try other three dynamical state, and check whether the parameter `g` and `ext_ratio`($\\nu_{ext}/\\nu_{th}$) you got consists with the phase diagram in Brunel's paper.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "execution": {}
   },
   "source": [
    "## End of tutorial"
   ]
  }
 ],
 "metadata": {
  "colab": {
   "collapsed_sections": [],
   "include_colab_link": true,
   "name": "W2_Tutorial2",
   "toc_visible": true
  },
  "interpreter": {
   "hash": "6c4039b0446a2e45c3d14928023f75989e3ca1f00bf719d898ca47d54d55585f"
  },
  "kernel": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "kernelspec": {
   "display_name": "Python 3.8.11 64-bit ('brian2': conda)",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.12"
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "window_display": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
